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1. INTRODUCTION 


Geomagnetic secular change has long been attributed to an imbalance 
between the effects of motional induction and magnetic flux diffusion 
within Earth's electrically conducting liquid outer core. In Part IA 
(VOORHIES, 1992) attention was focused on the fluid motion near the top 
of the core by adopting the source-free mant le / frozen- f lux core 
(SFM/FFC) magnetic earth model (wherein a rigid, impenetrable, 
electrically insulating mantle of uniform magnetic permeability 
surrounds a spherical, inviscid, perfectly conducting outer core in 
anelastic flow) . Several reasons were given to further consider the 
geomagnetic effects of motional induction by steady flow near the top of 
the core (e.g., a steady flow explicates quantitatively most of the 
recent, observed geomagnetic secular change). The theory underlying 
some estimates of core surface flow was summarized. Consequences of a 
few kinematic and dynamic hypotheses were derived: fluid downwelling is 
required to change the mean square radial magnetic flux density averaged 
over the surface of a FFC ; downwelling implies poleward flow for 
surficially geostrophic core motions. The solution of the forward 
steady motional induction problem at the top of a FFC was derived and 
found to be a fine example of deterministic chaos. Implications of 
persistent, if not steady, surficially geostrophic flow were described 
which apparently help explain certain features of the present broad- 
scale magnetic field and perhaps paleomagnet ic secular change. 

To investigate steady induction effects in geomagnetism, it is 
useful to regard the SFM/FFC model as a first approximation and to treat 
the supposition of steady surficial core flow as a hypothesis. To test 
hypotheses against observations, it seems appropriate to (a) understand 
both; (b) develop a satisfactory method for modeling the relevant 
observations in accord with the hypotheses; (c) apply the method to make 
quantitative predictions; and (d) subtract predicted from observed 
values and measure such residuals in units of the estimated uncertainty 
in the observations. In the context of the SFM/FFC approximation, this 
paper develops a method to fit the secular change indicated by 
geomagnetic field models in accord with the hypothesis of piecewise, 
statistically steady flow. Such field models represent the relevant 
geomagnetic observations very well and are here preferred to raw data 
for reasons noted in Part IA. 

With enough perfectly accurate information on the normal component 
of the time-varying magnetic flux density at the surface of a non- 
diffusive core, the steady surficial core flow inducing those variations 
could be uniquely determined by simple linear methods-provided the 
interior of the cryptic set is indeed empty (VOORHIES & BACKUS, 1985). 
It can be assumed that geomagnetic observations are not perfectly 
accurate and are sparsely distributed in space and time; therefore, they 
are not complete in either the spatial or temporal domain. It can 
further be assumed that models of such data are also imperfect and 
incomplete (e.g., truncated spherical harmonic models are incomplete in 
the spectral domain despite their completeness in the spatial domain) . 
It follows that there is not enough perfectly accurate information at 
Earth's surface, much less at the top of the core, to sustain the simple 
linear methods. One may, however, seek a steady surficial core flow 
which tracks that part of the total secular change indicated by real 
models of real data. 
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The inverse problem of deriving steady surficial core motions which 
fit imperfect models of observed geomagnetic field evolution is non 
linear due to the appearance of a transcendental exponential operator in 
the solution of the forward steady motional induction problem (IA) - 
Previous studies of the steady motions hypothesis (e.g., BLOXHAM, 
1987a, b, 1988a, b, 1989; VOORHIES, 1986a,b, 1987a, b; WHALER & CLARKE, 
1988) have noted the non-linear nature of the inverse steady motional 
induction problem. Methods for solving the non— linear inverse problem 
were developed and applied by VOORHIES (1987b) and BLOXHAM (1987b). 
Both methods feature iterative minimization of an objective function 
composed of: a square weighted residual relative to the secular change 
indicated by the geomagnetic field models fitted; an optional constraint 
requiring the flow to be as surficially geostrophic as desired; and an 
optional damping requiring the flow to be as smooth as desired. My 
method has, however, been based on a different approach which leads to 
differences in method, application, results, and interpretation. 

2. APPROACH 

For steady flow, the radial component of the induction equation at 
the top of a FFC of radius b is, in spherical polar coordinates (r,0,<}>), 

3 t B rp (b,t) + v s (b)«V s B rp (b,t) = B rp (b,t)d r u(b) . (0) 

This special case of the ROBERTS & SCOTT (1965) equation is but (13b) of 
Part IA. The subscripted p on the radial component of the magnetic flux 
density stresses that B rp is a prediction of the FFC approximation and 
the supposition of steady v s (b) during some time interval t 0 £ t £ t £ ; 
the extension to piecewise steady flow is straightforward. At the top 
of the free-streaming core (r = b, lb! =3,48 Mm) the components of the 
steady surficial fluid velocity v s (b) are still {u (b) =0 , v(b) , w(b) ] and 
V s « is still the surface divergence operator. Upward continuation from 
the core-mantle boundary (CMB) to Earth's surface (of B^ p (b,t) to 
B (a, t ) where lal = 6.3712 Mm) is also straightforward in the SFM 
approximation; however, an initial geomagnetic condition during [t 0# tf] 
is needed to use evolution equation (0). 

If the SFM/ FFC earth model were exact, and if complete and perfect 
knowledge of the radial geomagnetic flux density component at Earth's 
surface were available during some interval, then supposition of steady 
surficial core flow would overdetermine the inverse motional induction 
problem under otherwise fairly general circumstances (VOORHIES & BACKUS, 
1985) . A least-squares approach to solving this problem would then be 
formally justified and any residual misfit would falsify the 
supposition. In fact the SFM/FFC model is at best an approximation; 
moreover, the geomagnetic field is but imperfectly known at 'points' in 
space and time— so neither geomagnetic data nor spherical harmonic models 
thereof provide either complete or perfect information on the true 
radial field component. In view of the approximate nature of both the 
underlying physical assumptions and the models of limited data, I 
created and applied a weighted, optionally constrained, and optionally 
damped iterative weighted least squares method for solving the non- 
linear inverse motional induction problem posed by the hypothesis of 
steady flow. 
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The quantity to be minimized is the sum of three parts; a square 
weighted residual 4rcA r 2 measuring the misfit to the expected or input 
radial magnetic flux density at Earth's surface B r (a,t); an optional 
constraint 47eX g A 2 requiring the derived steady flow to be as nearly 
surficially geostrophic as desired; and an optional damping 4 XX 3 A 3 
requiring the spatial structure of the derived flow to be as smooth, or 
rather, as simple as desired. The total objective function is thus 

4nA 2 = 4n(A r 2 + XgAg 2 + *. d A d 2 ) (la) 

where X g and X^ are positive damping parameters. The magnitudes of \ 
and X^ determine respectively the importance of surficial geostrophy ana 
flow simplicity relative to the quality of fit. Large X g ensures 
surficial geostrophy; large X^ ensures simple flow. Two different types 
of weights are investigated: radial field weights and general weights. 

For radial field weights, the total square weighted residual 
accumulated during the interval from initial time t Q to final time tf is 
47 tA r 2 , where semi -normalized 


[A r (a;t 0 ,t f )] 2 s | <[B r (a,t) - B rp (a, t) ] 2 W(a, t)>dt 


(lb) 


<q(r,t)> denotes the mean value of q(r,t) averaged over the sphere of 
radius r (IA, equation (7)), and W(a,t) is the weight function. 
Clearly, A r 2 measures how poorly the predicted radial field B r p fits the 
expected (or input) radial field B r in units of the expected uncertainty 
OB r (a,t) a [W (a, t ) ] ~*/ 2 . For simplicity, the initial geomagnetic 
condition is taken to be B rp (a,t 0 ) = B r (a,t 0 ) with the understanding 
that a model for B r (a,t Q ) must not only be downwardly continued, but 
must be completed, to generate B r p(b,t Q ). If no SV were predicted, then 
B r p(a,t) = B r p(a,t Q ) = B r (a,t 0 ); then 4xA r 2 would be the total square 
weighted change of the radial field accumulated during the interval (the 
square weighted signal). If SV were also constant, then 4xA r 2 / I tf-t c I 
would increase only if the weight function grew heavier with time; 
therefore, A r 2 is normalized for sphericity but not interval. More 
generally, the fit would be judged adequate if A r 2 £ ( e *U*/ if 

The integrand in (lb) is the 


W ( a , t ) “ 1 ) 


[B r (a,t) - B rp (a,t ) ] 2 ... 

instantaneous square weighted residual [ 8 r (a, t ;t Q ) ] 2 
B rp (a,t) ] 2 W(a,t)>. 

The mean square ageostrophy of the flow 


<(B r (a,t) 


[A g (b; t Q , t f ) ] 2 s <[d r u(b)cos0 + v(b) sin0/b] 2 > (lc) 

measures departures from the geostrophic radial vorticity balance (IA, 
equation (12) whereby downwelling (3 r u > 0) implies poleward flow). The 
geostrophic radial vorticity constraint is not needed to derive formally 
unique, piecewise steady core surface motions, but it is plausible 
dynamically. In the limit as X g approaches infinity, this constraint is 
consistent with tangential geostrophy-which eliminates the toroidal 
ambiguity in B r v s (BACKUS, 1982) in some areas (BACKUS & LEMOUEL, 1986; 
HILLS, 1979) and reduces it everywhere on the CMB. More generally, this 
constraint reduces the geomagnetic information required to uniquely 
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determine a steady flow. Surficially geostrophic flows have also been 
used to estimate the purely mechanical or topographic torque exerted by 
the core on the mantle (SPEITH et al . , 1986) and, with the added 
supposition of tangential geostrophy, the perturbation pressure field at 
the CMB (VOORHIES, 1991) . Occasional application of this interesting 
constraint needs no further justification. 

With radial vorticity tOj. = r»Vxv and anelastic downwelling 9 r u(b) = 
-V *v (b) , the measure of spatial complexity adopted is the sum of the 
mean square radial vorticity and the mean square downwelling of the flow 

[A d {b;t 0 ,t f ) ] 2 ^ < [d) r (b) ] 2 + [3 r u(b)] 2 > . (Id) 

By ( lb-d) , both X,g and X d in (la) must have dimensions of time cubed. 
Finitude of (Id) ensures piecewise continuous fluid velocity v s (b) ; 
truncated spherical harmonic models of v s (b) are continuous and smooth. 

A plausible (sub-relativistic and sub-acoustic) core surface flow 
need not be spatially simple, nor are very smooth flows necessarily more 
reliable than other flows. The bias towards simple flow was introduced 
chiefly to speed convergence of the iteration scheme. It turns out that 
varying X d enables exploration of how well various steady flows fit 
secular change. For example, if the SFM/FFC earth model and steady flow 
admit A 2 5 I t f -t Q l , then A. d may be chosen so as to achieve an adequate 
fit and eliminate unnecessary spatial structure in the flow; if not, 
varying X d allows a reckoning of how much SV might reasonably be 
attributed to a steady flow at the top of a FFC surrounded by a SFM. 
Indeed, X d might be chosen to give the simplest flow yielding the 7% 
residuals expected in the SFM/FFC approximation (IA). Such a choice 
might reduce diffusive effects on the estimated downwelling anticipated 
by MUTH (1984, pers . comm.) and BLOXHAM (1989). Moreover, non-zero X d 
requires both the radial vorticity and the downwelling to be finite when 
averaged over any non-zero area-in partial accord with the condition of 
hydrodynamic motion (BONDI & GOLD, 1950) . Non-zero >. d further ensures a 
surface kinetic energy density spectrum (VOORHIES, 1986a) which falls 
off faster than n -2 at very high spherical harmonic degree n. Then the 
total surface kinetic energy density will converge with n— even in the 
absence of geomagnetic information resolving small-scale flow structure. 

Other constraints could be substituted for or added to A d . One 
tempting constraint imposes both finite downwelling and finite radial 
vorticity at all points b. This would prohibit vortex sheets, vortex 
lines, and singular downwellings just within the fluid core-in full 
accord with the condition of hydrodynamic motion. Any core analogs of 
fronts, tornadoes, plumes, or boundary currents would then have finite 
thickness-even lacking the geomagnetic information needed to resolve 
them. Yet tearing of the fluid is not implied by the use of truncated 
spherical harmonic representations of a finite fluid velocity field, nor 
does estimation of a truncated parameter set imply that unestimated, 
higher degree coefficients are zero. One might place upper bounds on 
fine structure thickness and reduce ringing by extending the maximum 
degree of such an estimate far into the damping regime (wherein damping 
rather than geomagnetic information and non-zero molecular dif fusivities 
determines the scale of fine structures) . This may be too burdensome 
computationally. Moreover, as a norm (Id) is consistent with 
effectively inviscid flow; norms barring sheet vorticies seem 
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inconsistent with (0). So (Id) offers a good compromise between 
smearing out any fairly sharp-edged jets, gyres, and plumes needed to 
fit secular change and ringing caused by truncation. Furthermore, if 
localized diffusive flux eruption or decay masquerade as strong frozen- 
flux upwelling or downwelling plumes, then the milder constraint (Id) 
allows, and indeed encourages, spatial confinement of such artifacts. 

Arbitrary selection of X^ * 0 injects prior bias rather than 
genuine prior information into what some might otherwise view as a 
Gauss-Markov estimation problem. Following BACKUS (1988a), the damped 
weighted least-squares approach is not stochastic inversion nor is it 
properly Bayesian inference when X^ is varied to investigate various 
flows rather than impose, a priori, a particular personal probability 
distribution upon the flow parameters. Unfortunately, with arbitrary X ^ 
* 0 the derivation of reliable uncertainty estimates for the velocity 
field parameters is difficult or impossible. If the contribution from 
X^A^ to the total information matrix were replaced by genuine prior 
information before inversion, then the resulting covariance could be 
physically meaningful. Prior information on core motions includes: (i) 
the time-averaged viscous dissipation within the core must not exceed 
the time-averaged geothermal flux; and (ii) the core flow speeds of 
interest must be less than Mach one everywhere and are likely much less 
than the rotational speed (253 m/s at b = 3.48 Mm and 0 = n/ 2). The 
former places no constraint on effectively inviscid motions at the top 
of the core; the latter is too weak to speed convergence of the 
iteration scheme. Yet the formal uniqueness problem is solved, the 
existence problem is of immediate interest, and the question of 
practical uniqueness within reliable uncertainty estimates is moot if 
existence cannot be established. When seeking plausible solutions to 
the existence problem, baseless bias towards smooth flow ought not 
hinder hypothesis testing. Such bias can be reduced (or eliminated) by 
reducing X ^ towards (or to) zero or by modifying the algorithm as 
described in section 3.3. 

An alternative form for the square weighted residual (lb) suitable 
for use with 'discrete' weighted geomagnetic data D(rj,tj) (be it D, I, 
H, X, Y, Z, or F) gathered between radii r^ and r Q is 


t f r Q 2n it 

d I / J [D(r,t)-D p (r,t)]jM.j k [D(r,t)-D p (r,t)] k r 2 sinededif>drdt 


to r i° 0 


where Mj^ is the appropriate weight matrix function reflecting any 
expected correlation between the j th and k^* 1 data and C is the 

■j *5 

appropriate semi -normalization constant {3/47C[r Q -r^ ] } (Voorhies, 1988 
unpublished) . LANGEL (1990, personal communication) stresses that this 
form leads to simultaneous estimation of both an initial geomagnetic 
field model and a steady surficial core flow. Although the formalism 
and preliminary solutions lie outside the range of the present series, 
it can be seen that the combination of secular change data with the 
steady motions hypothesis places powerful, albeit perhaps contrived, 
constraints on an initial geomagnetic field model at the CMB. 
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3. METHOD 


3.1 The Square Weighted Residual 

With initial condition B r (a,t 0 ) = B rp (a,t 0 ), (lb) is rewritten as 


t f t 2 

A r 2 = J <{J [8 x B r (a,x) - a t B rp (a,T) ]dT) W(a,t)>dt (2) 

fc o fc o 

where the dependence of A r 2 on (a;t Q ,tf) is understood. For r ^ b, the 
input radial field B r (r,t), the predicted radial field B rp (r,t), and 
their time derivatives are expressed in terms of their compact spherical 
harmonic expansions 

B r (r,t) = g i (r,t)S i (0,^) 3 t B r (r,t) = 3,^ (r, t) Si (0, + ) (3a, b) 

B rp (r,t) = Y i (r,t)S i (0,^>) 3 t B rp (r,t) = 3 t Y^ (r, t ) (0, <j>) . (4a, b) 

Repeated subscripts are summed over (Einstein convention); the spherical 
harmonics S ± ( 0 , <(» ) and radial field coefficients are defined as follows. 
Let P *** represent the Schmidt normalized associated Legendre function of 
degree n and order m (CHAPMAN & BARTELS, 1940; JACOBS, 1987). For index 
i = n 2 , S i (0,<^) « P n °(cos0); for i = n 2 + 2m-l and m*0, s i( 0 , <t> ) = 

[cosm<fr]P n m (cos0) ; for i = n 2 + 2m and m*0, S^O.^) = (sinm<|>] P r 7 n (cos0) . 

Clearly, n(i) and m(i) are specified by i. In the SFM approximation, 
the expected radial field coefficients g^r.t) are the corresponding 
input Gauss coefficients (g n m , h n n ') multiplied by [n+1] [a/r] . The 

predicted radial field coefficients y i (r,t) are similarly defined and 
are derived via spherical harmonic analysis of the radial field 
predicted by steady motional induction at the CMB for t ^ t Q . Radii a 
and b are of primary interest, so let 


g i (a,t) = g^t) = 
Y i (a,t) = Yi<t> = Yi 

«i = Y ij G i 


gi (b,t) = G^t) = G i 
Yi (b, t) = ^(t) = Ti 

^ " Y i j F j 


(5a, b) 
(6a, b) 
(7a, b) 


w here the time dependence is understood and, with 5^j denoting the 
Kroenecker delta, the upward continuation operator is a diagonal matrix 
with elements = [b/a] n ( +2 8i j . For a SFM. the map the radial 
field coefficients, and thus the scaloidal core field, from the CMB to 
Earth's surface. Curiously, a diagonal upward continuation filter with 
elements that depend upon n and the temporal frequency of the 
discrete Fourier transformed g^(r,t) -4 g^ (r , (0^) can account for the 
effects of non-zero, laterally homogeneous, mantle electrical 
conductivity; however, mere inclusion of m-dependent and off-diagonal 
elements will not account for the toroidal-poloidal coupling expected 
for laterally heterogeneous mantle conductivity (Voorhies, 1988 
unpublished manuscript) . 

The s t ream func t ion -T(b) and the velocity potential -U(b) 
describing the steady surficial fluid velocity field v(b;t 0 ,tf) (hence 
B r p(b,t*t c )) are also expanded in terms of spherical harmonics: 
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▼ (b) = V g T(b, 0,<(>)xr + V s U(b,e,<(» (8a) 

T (b) = a i S i (Q t $) U (b) = 0^(0, $) . (8b, c) 

This ensures V*v s = 0 at b-as required by the kinematic boundary and 
anelastic flow conditions (IA). 

Spherical harmonics are orthogonal on spheres, so (3a-7b) can be 
used to rewrite (2) as 



t f 

<[ g i s i 


2 

YjSi] W(a,t)>dt 


(9a) 


= J <(9i ~ Yj.) [SiW(a,t)Sj] (gj - Yj)>dt 


or 


4jcA r 2 = J [ gi - Yj. 1 w i j (9j ' 7j]dt 

t„ 


(9b) 


(9c) 


where the dependence of the on (0,$) is understood and W^j(t) s 
4rc<S^W(a, t) S j> defines the elements of the time-dependent radial field 
weight matrix. In the case of equal weighting (VOORHIES, 1986b), W— 
reduces to the diagonal normalization matrix for Schmidt normalized 
spherical harmonics, N^j ^ 4n [2n ( i ) +1 ] j . If the g^ describing 
B r (a,t) are not all equally well determined, not independent, or both, 
then W^j does not reduce to 

For general weights, // replaces W^j/4JC in (9c): with matrix 
inversion preceding assignation of element indices, * 
Eo defines the time-varying covariance associated with the input 
radial field model; the e^ are the unknown true errors in the input g^; 
and the Eo operator yields the expected value (see Appendix) . 

Radial field weights rely upon a scalar weight function W(a,t) 
which is the inverse of the expected squared uncertainty in B r (a,t). If 
uncorrelated observations of B r were used, then W would be the inverse 
squared uncertainty of the observations when and where observations 
exist; W would be zero elsewhere and elsewhen as non-existent data enjoy 
zero weight. When spherical harmonic models of the radial field are 
used the weight function is 

W(a,t) s [OB r (a, 0, <|>, t) ] -2 = (S k (e,<j»)E kl (a,t)S 1 (0 / 4O] _:l . (10) 


If the covariance for the Gauss coefficients at (a,t) is V k j, then 
E kl (a,t) = [n(k)+l]V kl [n(l)+l] » R ik V iH R ilr where R ik * (n(i)+l]o ik . 
Substitution of (10) into (9b) yields 


A r 2 = J f <( 9i - VV^EjaSii’^gj 

fc o 

or 


Yj ) >dt 


(11a) 
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(lib) 


tf 

4nA r 2 = J I (g ± - Yi)Wij(gj - Yj)]dt 
to 

tf t t 

= 1 { [J O x 9i - 3 t Yi)dX]W i;j [J O t g-i - 3 x Yj)dX])dt 

to to fc o 

where the radial field weight matrix elements are 


-1 

W i j(a,t) = 4re<S i [Sj^Ej^SjJ Sj> 
= 4rt<S i [w k (a, t)S k ]Sj> 

= 4 jtZ ijk w k ' 


(12a) 

( 12 b) 

(12c) 


the time-dependent spherical harmonic representation of the weight 
function (including its non-zero mean value) is 


W(a,t) s w k (a, t)S k (0,4>) 


[S k E k i( a ,t)Si ]- 1 


(13a) 


and the symmetric third-rank tensor has elements 

z ijk s <s i s j s k > = z jik = z ikj ■ (13b) 

When the weight function is independent of position, then is 

proportional to . If W(a,t) is everywhere and always equal to unity, 
then W- • reduces to N ii . The latter conditions were in effect presumed 
by VOORHIES (1986b); such presumptions are avoided here. Although 
W(a, t) should be nearly laterally homogeneous for broad-scale models of 
satellite data, such data are not always available; when they are, 
W(a,t) can be (16 nT )' 2 (LANGEL, ESTES, & SABAKA, 1988a; 1989). A 

derivation of weight matrices for the Definitive Geomagnetic Reference 
Field (DGRF) models ( IAGA, 1988) will be described in Part IC, 

If SV coefficients were used as input and if an expected error 
covariance for the time rate of change of the radial field coefficients 
(E®^ a Eot (df-ej (d,.e H )}) were available, then (11c) would become 

4rtA 2 = r f ur 4i<a T gi - d x Yi>*nJ t 4jO X 9j - d X Yj)<at])dt 
c o c o 0 

where the time (T) -dependent matrix © SV is the upper triangular m|trix 
square root of the SV weight matrix defined as in (12a) but with E k ^ 
replacing E kl . This approach was not pursued. Though the DGRF models 
employed were used to derive dummy SV models for the x integration, main 
field weights seem more appropriate to the fitting of a sequence of main 
field models and were thus used to weight the residuals. 

In order to minimize the objective function with respect to the 
flow parameters 04 and P i( A r 2 in (He) must be expressed in terms of 
these parameters. By analogy with VOORHIES (1986b), write (0) as 

3 t B rp (b,t) = -V s *[B rp (b,t)v g (b)] < 14a) 
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(14b) 


or, using equations ( 2 ) through ( 8 ), 


3 t r k s k 


= -V< r i s if V 3< a j s j) xr + 


= b-^t 


d 4» s i a 0 s j 

sin 0 


3 0 s i a (|> s j 

sin 0 




+ 


b _ 2 r" i {Si [n ( j ) ] [n ( j ) + 1 ] Sj - deS^Sj 


ViVj 

sin 2 0 


>Pj 


(14c) 


Left multiply the scalar equation (14c) by S^sin0/4rc and integrate over 
0 and <|> . Then use the orthogonality of the spherical harmonics to 
evaluate the left-hand side, relabel 1 -»k, and reorder the integrand on 
the right-hand side noting terms like S^SjPj = SjSjJJj because there is 
no sum over k. The result is written 


^t^k 58 ^i x ijk^ a j + ^i Y ijk^Pj " p kj a j + ^kjPj 

where to develop the third-order coupling, the third-rank tensors are 


x ijk 

Y ijk 


b' 2 [2n(k)+l]<O ( | ) S i a 0 Sj - aQS^Sj }S k /sin0> = 

b - 2 [2n(k)+l]<{S i [n(j) ] [n ( j ) +1] Sj - deS^Sj 


~ X jik 


ViVj 

sin 2 0 


)S k > 


(16a) 

(16b) 


and, 

p kj 

Q kj 


for numerical integration over the CMB, the second-rank tensors are 


b" 2 [ 2 n(k)+l]<S k {a ( j ) B rp d 0 Sj - deBj-pd^Sj } /sin 0 > 


(17a) 


b -2 (2n (k) +1] <S k (B rp [n ( j ) ] [n ( j ) +1] S j - deB^Sj 


VrpVj, 

}> • 

sin 2 0 (17b) 


Upward continuation of (15) via (7b) yields 
d t 7i = ^*ik^t^k = ^ik^^i x ijk^ a j + ^i Y ijk^Pj^ 

= Yi^lPkjOtj + Q^jPj] 

5 Pij a j + qijpj - A n ^ . 

In (18c) Pij a ' r ik p kj ; ^ij s T ik Q kj ; ^1 * (Oj.’Pj) defines the extended 
vector obtained by concatenating the coefficients of T(b) with those of 
U (b) ; and time-varying A^ = ^Pij ;c *ij} defines the extended matrix 
obtained by concatenating the corresponding sub-matrices. Substitution 
of (18c) into ( 11 c) and relabeling yields 

t f t t 

4nA r 2 = J ([J 0 T gi ~ A ik 5 k )dt]W i;j [f - A jk 5 k )dT] )dt . (19) 


(18a) 

(18b) 

(18c) 
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The square weighted residual (19) will be minimal only if it is 
extreme, in which case its total derivatives with respect to the 
parameters ^ vanish. The first approximation to this condition was to 
set the partial derivatives of (19) with respect to zero. Then 


t f t 


Oj = -2 S ([/ < A ik 5 kl )dT]w ij (3 x g j " A jk^k )dtl)dt 


(20a) 


fc o ^ 


where the ultimate dependence of the A^ upon B r p(b,t), hence on v(b) 
and the initial condition (i.e., on and G^(t Q )), has been temporarily 
omitted to achieve linearization. (Recent tests show this is a good 
approximation, particularly for slow flows over short intervals) . In 
matrix notation equation (20a) is 


tf t t 

0 = -2j { (J A dt] T W [J (3 x g - A £) dx ] } dt 

t Q t Q = “ t Q 

where a single underline denotes a column vector, a double underline 
denotes a matrix, and the T superscript indicates the transpose. The 
linear least-squares (LLS) estimate of the parameters is 


t f t t -i 

^LLS _ (J ( [J A dX] T W [J A dx ] } dt ) 

to fc o = = fc o = 

t f t t 

{/ { [j A dx] W [J 3 T g dx] }dt) (20b) 


= (A T W A) _1 (A T W 3 x g) (20c) 

where the lower (single) overbar indicates dummy time integration over X 
from t to t, and the upper (double) overbar indicates total time 
integration over t from t Q to t^. Equation (20c) provides the linear 
least -squares estimate of the steady streamfunction and velocity 
potential coefficients which best fit the expected evolution of the 
radial magnetic flux density at Earth's surface during the interval. 

There are three crucial differences between (14a-20c) above and 
equations (10-17) of VOORHIES (1986b). Firstly, the squared residual is 
non -uniformly weighted. Secondly, the square weighted residual includes 
the double time integral needed to fit the evolving main field rather 
than simple SV. Thirdly, the predicted radial field (B rp or its 
coefficients T i ) appears on the right-hand sides instead of the input 
field (B r or its coefficients G^) . The A^^ in (19) thus depend upon the 
predicted main field at the CMB instead of the input field model. The 
elements of A in (20) also depend upon the predicted field at the CMB 
instead of tHe input field, but this dependence was suppressed to obtain 
a system of linear equations for the model parameters |. 
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To solve the non-linear problem, I developed an iterative method 
wherein the linearized problem is first posed by supposing for 

purposes of estimating the elements of A. These can be computed using 
equations (15), either (16) or (17), and (18), The linearized least- 
squares problem (20c) is then solved for £ and thus v(b) by back 
substitution into (8). B r p(b,t) is computed by solving the forward 
motional induction problem (14a) from the initial condition B rp (b,t 0 ) = 
B r (b,t Q ). These values for B r p(b,t) are then used to compute r and 
residuals (g - *)[) ; both weighted and unweighted residuals; and new 
elements for A matrix elements via (15), either (16) or (17), and (18). 
The new A matrices and the residuals comprise the input for the next 
iteration. Let j indicate the such deep matrix iteration, let 
|(j+l) = £(j) + 8§(j + l), and let 5g(j + l) = g - ]f(j) so that f&^gfj + l) = 
(d x g_ - 3 x y( j ) ) • Then 


gj:NLLS ( j + X ) 


(A T (j) W A( j ) ] 


-1 


[A T ( j ) W 5d x g(j) ] 


( 21 ) 


describes this non-linear least-squares iteration procedure. 

For t*t a the spherical harmonic content of B rp typically extends to 
far higher degree than that of B r ; such narrow-scale structure must be 
preserved in the computation of the new A, though it will typically not 
contribute to the weighted residuals due”to truncation of g^, hence j . 


3.2 The Geostrophic Radial Vorticity Constraint 

Parameterization of the geostrophic radial vorticity constraint 
proceeds by writing the zero-mean geostrophic radial vorticity balance 
at the top of the core in terms (e^) of its spherical harmonic expansion 


v 

e k s k = (3 r ucos0 + -sinB) 


(22a) 


= b^tOtj^Sj + (3 j {3 qS jsin0 + n ( j ) [n ( j ) +1] Sjcos0) ] . (22b) 

Left multiply (22b) by S^sin0, integrate over a spherical CMB, and use 
the orthogonality of the spherical harmonics to obtain 

c i = C ij a j + D ij^j * B ik^k (23a) 

where, corresponding to is the appropriate element of either 

C i j * b~ 2 [2n(i)+l]<S i 9 ( | ) Sj> (23b) 

or 

D^j s b“ 2 [2n(i) +1 J <S^ {9gS j sin0 + n ( j ) [n { j ) +1] Sjcos0}> . (23c) 


Only spherical harmonics of like degree and order contribute to the 
j . Indeed, only certain elements adjacent to the diagonal of matrix C 
(corresponding to cosm<|> O^sinm^) ] and sinm<|> [3^cosm4>] ) are non- trivial; 
these are readily evaluated analytically. Only spherical harmonics of 
like order contribute to the D^j . However, if the spherical harmonic 
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expansions of streamf unct ion -T(b) and velocity potential -U(b) are 
truncated at degree and order Nrp = Ny = N v , then ^ij (hence B.^) can 
have non-trivial elements for i < [N v +1] [N v +3] - whicn can include i > 

N [N +2] . This is due immediately to the factors of sin© and cosQ 
appearing in (22) and ultimately to the latitudinal variation of the 
radial planetary vorticity = 2£^cos0 . Neglect of this fact (e.g., 
assuming = 0 for i > N v (N v +2)) may lead to a truncated velocity 

field which fails to be even approximately surficially geostrophic. 

Equations (22) and (23) allow the geostrophic radial vorticity 
constraint (lc) to be rewritten 

4irXgA g (b; t Q , t f ) 2 = 47t>.g<lE i s i ] Ce-jS-j )> 

= ^gf e i N ij E j^ 

= (B ^ ) ^g ( B 5) 

where \ Q is the normalization matrix multiplied by X g . Adding (24d) to 
(20c) gives the constrained, weighted objective function 


47t[A r 2 + X g A g 2 ] = (d x g “ A ^) T W(3 T g - A £) + ^) T A g (B . (25) 

This function is minimal only if it is extreme, in which case its 
partial derivative with respect to any element of £ vanishes. Again 
omitting the weak dependence of A on 2^, the constrained, weighted linear 
least-squares estimate of the parameters is 


(24a) 

(24b) 

(24c) 

(24d) 


jiCLLS u {A T w a + b T A g B) 1 (A T W 3 t g) . 


(26) 


In the context of the iterative solution to the non-linear, 
constrained, weighted inverse problem 


g^CNLLS ( j + 1) _ [ A T ( j) w A( j) + b T A g B] 1 {A T (j) W 5d T g ( j ) 

- [f A g B][*(j) - 5g 0 ]} (27) 

replaces equation (21) . The prior estimate of the parameters is | go . 

If £ =0 the bias is fixed on zero mean with intensity Ag. If |g 0 = 

|(jj^ 9 the bias floats with the iteration scheme-yielding a learning 
algorithm. In the latter case, small X g may yield large corrections 
8£(j+l) which result in a rather ageostrophic velocity field far from 
that used to calculate A(j) on iteration j; one might prefer an 
iteration-dependent constraint parameter with such a learning algorithm. 
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Here attention is focused on J; go = 0 and values of X g which are 
either zero for the unconstrained piecewise steady inverse motional 
induction problem or so large that the normalized mean square 
ageostrophy (VOORHIES, 1986c) 

_ < [3 r ucos9 + vsin0/b] 2 > 

9 = <[a r ucos0] 2 + [vsin0/b] 2 > 

is less than 10~^. In the latter case, the constrained weighted least- 
squares solution will equal that obtained using stochastic inversion 
with prior information matrix B T A g B. Because the geostrophic constraint 
is viewed as a plausible approximation rather than an absolute necessity 
required by genuine prior information, solutions (26) or (27) are not 
considered to be stochastic estimates as described by MCLEOD (1986). 
These estimates involve 'soft bias' rather than soft or hard bounds as 
described by BACKUS (1988a,b) . A hard version of the geostrophic radial 
vorticity constraint might be imposed using the geostrophic basis 
functions of BACKUS & LEMOUEL (1986). 

Intermediate values of X g can be used to study the tradeoff between 
the constraint and the length of the interval in which the flow is 
assumed to be steady. This seems appropriate when the surficially 
geostrophic flow hypothesis is tested in the context of a SFM/FFC model 
with piecewise steady surficial core motions. This approach was taken 
using the unweighted, non-iterative linearized method of VOORHIES 
(1986b). Analysis of the resulting, severely truncated (N v = 5) fluid 
velocity fields (VOORHIES, 1986c) showed that the constraint reduces the 
tightness of fit; yet a mild constraint can increase slightly the 
accuracy of geomagnetic forecasts made by extrapolating the effects of 
steady motional induction at the CMB outside the interval t Q < t < tf 
and subsequent upward continuation. The former point has been confirmed 
using superior methods [VOORHIES, 1987c, 1988, 1989; BLOXHAM, 1988b]. 

3.3 Damping Mean Square Radial Vorticity and Mean Square Downwelling 

If the SFM/FFC earth model and the working hypothesis of steady 
(optionally geostrophic) flow were correct, and if complete, albeit 
imperfect, information on the evolving geomagnetic field were available, 
then the weighted (constrained) least-squares estimate (26-27) would 
uniquely determine the (constrained) steady flow to within uncertainties 
implied by the inverse of the total information matrix. This covariance 
of the flow parameters [A T WA + B T A g B] is mere expectation; it depends 
upon neither the flow parameters” nor the residuals. Because complete 
information is not available, if one steady flow adequately fitted the 
incomplete, imperfect information, there may well be others which do so. 
In seeking whether one such flow exists, it seems reasonable to 
initially restrict attention to solutions which are spatially simple. 

I chose to seek flows characterized by low values of the mean 
square surficial curvature of both the streamfunct ion and the velocity 
potential in hopes of eliminating unecessary flow structure- 
particularly on small spatial scales. The radial vorticity 

A 

<o r (b) = r*Vxv g (b) = -V s 2 T(b) = a^lnd) [n(i)+l]S i ) (28) 

is the surface Laplacian of the streamfunct ion; the downwelling 
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(29) 


a r u(b) = -Vg*v 3 ■ -Vg 2 U(b) * P^nd) [nUl+USjJ 

is the surface Laplacian of the velocity potential. A velocity field 
with small mean square radial vorticity and small mean square surface 
divergence has small mean square surface curvature in both T and U. 
This choice will tend to fill in any regions where the velocity field is 
relatively poorly determined by interpolation from surrounding regions 
without introducing unecessary flow sources and without smearing out 
isolated jets, gyres, or plumes. 

Parameterization of equation (Id) using (28) and (29) yields 

4jlX d A d (b;t 0 ,t f ) 2 = 47& d <[V s 2 T] 2 + [V g 2 U] 2 > 

~ a i F ij a j + ^i F ij^j * \ A d \ 

where F Aj * 8^^ j4nX d (n ( i ) ] 2 [n ( i) +1 ) 2 / [2n ( i) +1] and is the extended 
diagonal matrix with elements in both the upper left and lower right 
quarters equal to those of Fjj . Either <[(0 r (b)] 2 > or <[3 r u(b)] 2 > may 
be damped more strongly by adjusting the elements of This option 

was not pursued despite earlier findings (VOORHIES, 1584, 1986a) 
suggesting more energetic toroidal flow. 

Adding (30c) to (25) gives the damped, constrained, weighted 
objective function (1) 


4rt[A r 2 + X g A g 2 +X ( ^A c j 2 ] ^ (d^ A £) T W (d T g A ^) + 

(B 4) T A g (B \) + (f A d ^) . (31) 

This total objective function is minimal only if it is extreme, in which 
case its total derivatives with respect to the elements of % all vanish. 
The weighted, constrained, and damped linearized least-squares estimate 
of the parameters puts the partial derivatives to zero and is given by 


(30a) 

(30b) 


£ a (A T W a + B T Ag B + A d ) _1 (A T W 9 t g) 


(32) 


In the context of the iterative solution to the weighted, 
constrained, and damped non-linear inverse problem, 


8£( j+1) = [A T ( j ) W A(j) + B T Ag B + 
- [B T Ag B] [!j(j) - ^g Q ] 


A d 5_1{AT( 3) w 53 x g(j) 
- A d [^_(j) - s do n 


(33a) 


replaces equation (21) or (27). 
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The value of deep matrix element iteration depends upon how much 
A ( j ) changes with j. This in turn depends upon how incompatible the 
Input field models, hence A(j=0), are with the earth model, and upon the 
intitial condition, the estimated fluid velocity, and the length of the 
interval during which steady flow is presumed. Faster flows tend to 
generate appreciable non-linear feedback more rapidly; yet even a slow 
flow may do so eventually. Rough flows tend to generate small-scale 
structure in B r p(b,t) via a chaotic cascade to ever higher wavenumbers. 
This need not degrade the fit to broad-scale input field models because 
such models do not specify small-scale field structure (i.e., because 
unknown high-degree Gauss coefficients are assigned zero weight). In 
fact, high-degree structure in B r p(b,t) can be exploited by the flow to 
improve the fit to the broad-scale field models (via a reverse cascade) . 
Unlike broad-scale spherical harmonic models, surface data may be used 
to test (or constrain) high-degree structure in B (a,t); measures of 
B r p(b,t) and its time rate of change showed that the steady flows 
derived from the DGRFs were not rough enough to grossly violate such 
constraints during the several-decade interval targeted. 

As noted above, ^g Q was taken to be 0 in the actual calculations. 
Though is commonly taken to be 0, many calculations have been 

performed with floating bias — »^(j) and an adjustable convergence 

factor X^ — »^ d (j+l). The resulting learning algorithm is useful for 
deriving rougher flows which otherwise require values of fixed X^ so 
small as to inhibit convergence of the iteration scheme or even allow 
convergence towards local minima where A r ^ exceeds values found using 
larger X^ or the learning algorithm. The learning algorithm thus 
relaxes restrictions imposed by otherwise-f ixed bias toward zero flow 
roughness on tests of the steady flow hypothesis. Flows derived using 
the learning algorithm are, of course, not optimally simple; they can be 
smoothed by switching to the fixed X^ iteration scheme with bias towards 
zero flow roughness. Some may prefer this strategy to imposing both a 
fixed bias toward no flow and a floating bias toward the previous 
estimate governed by a convergence factor. 

The use of a non-trivial flow estimate £{j) to predict B r p(b,t) and 
calculate new A(j) may seem inconsistent with a fixed bias towards 
parameters "which are zero; however, the fixed bias strategy need not 
be inferior to the learning algorithm; indeed, the former may yield 
flows which represent any steady part of the true flow near the top of 
the core more accurately than those derived by excessive repetition of 
the learning algorithm due to errors in the SFM/FFC approximation. 

Because deep matrix element iteration can be computationally 
burdensome when equations (16) are integrated numerically over the CMB, 
it seems worth ensuring that the best estimate of £(j) is used to obtain 
B rp t) for the subsequent calculation of A( j + 1) . In principle, this 
can be accomplished by introducing shallow iteration whereby the 
correction vector of streamf unction and velocity potential coefficients 
determined on sub-iteration i+1 of deep matrix element iteration j is 


-1 


5^ (i + 1, j) = [ A T ( j ) W A ( j ) + B t Ag B + A d ) { A t ( j ) W 53 x g(i,j 
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(33b) 


- [B t A g B] [jj ( i, j ) - ^ 0 1 - A d [^_(i, j) - §do ]) • 


Shallow iteration proved to be of but slight use in practice. 

4. SUMMARY 

In the SFM/FFC earth model derivation of a (piecewise, statistically) 
steady fluid flow near the top of the core from imperfect models of the 
time-varying geomagnetic field requires solving a non-linear geophysical 
inverse problem. In 1987, a method was developed to solve this problem. 
The method attempts iterative minimization of an objective function 
composed of the square weighted residual in the geomagnetic secular 
change relative to a reference epoch; the ageostrophy of the flow as 
measured by the mean square departure from a geostrophic radial 
vorticity balance; and the spatial complexity of the flow as measured by 
flow source amplitudes (the mean square radial vorticity and the mean 
square downwelling of the flow) . The geostrophic constraint and the 
damping of spatial structure are optionally imposed with variable- 
strength damping parameters. In order to mitigate the artificial 
restrictions imposed by prior bias towards zero flow sources on the 
investigation of steady flows, a learning algorithm was also developed 
in which the bias is shifted towards the result of the previous 
iteration and departures therefrom governed by a convergence factor. 

When combined with numerical forward solution of the surficial 
motional induction equation (14a), equation (33a) defines a method for 
solving iteratively the non-linear geophysical inverse problem posed by 
the simple suppositions of a source-free mantle surrounding a frozen- 
flux core in surficially steady motion. The method involves two levels 
of iteration, double time integrals of matrices whose elements are 
surface integrals, a weight function which varies with both position and 
time, and includes two optional constraints: one for imposing the 
geostrophic radial vorticity balance and one for damping the spatial 
complexity of the flow. Additionally, two kinds of weights have been 
explored with fixed bias and learning algorithms. The reader may 
appreciate the practical difficulties of applying this method and 
keeping track of the various solutions, the irony of having such a 
simple set of working hypotheses yield so complicated a method, and most 
importantly, the complexity of the real Earth when stripped of such 
simplifying suppositions. 
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APPENDIX 


Consider the use of general weight matrix elements // instead 
of W- -/4 ji in (9c) . The generalized weight matrix £2 is the inverse of 
the error covariance matrix for the radial field coefficients E: 


hj a Eo ^ e i e j ) 


^ij E jk ” ^ik 


(Ala) 


(Alb) 


where e^ is the (unknown) true error associated with the use of expected 
radial field coefficients g.^ and only regular expectations (non-singular 
E) are considered, 
weighted residual 


This transforms A.. 2 into the generalized square 



(A2) 


. fc f T 

= J (g - y) £2 (g - Y)<*t 
t Q " = 

where matrix notation is employed in the last step. The scalar /( t) 
renormalizes the instantaneous generalized square weighted residual 
[8 *(a,t;t Q )] 2 , hence A r * 2 . It is taken to be the trace of fl(t)E(t) 

— typically the number of radial field coefficients fitted at time t. 

No correction for the number of degrees of freedom of the flow model is 
included because it is the significance of the residuals, rather than 
the efficiency of the model, which is of interest here (see Part IC) . 

Let co k j denote the elements of the upper triangular matrix square 
root of the positive definite generalized weight matrix 


ft = C0 T G) 


(A3) 


(BIERMAN, 1977, p40). Then 


/[A r *(a;t 0 ,t f )] 2 = 


t f T t 

J (g - y) to oo(g - Y)dt 


T T 


-1 


J 4jc< (g - Y) <0 (S N S) 
t„ " - = - - - 


a> (g - Y)> dt 


= J 


t f 

4n<[B r (a,t) 


2 

B rp (a,t)] >dt 
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(A4) 


, t f ~ - - 

=J (g k - Y k ) (g k - Y k > dt 

fc o 

where 

2n(k) +1 1/2 - - 2n(k) +1 1/2 

B r -g k ([ ] S k ) B rp S7k ([ ] s k ) 

4 n 4rc 

and where g k ® ©j^g.^ or y k s © ki y i are respectively the input or 
predicted radial field coefficients measured in units of expected 
uncertainty . 

To better understand the general weights, introduce the vector of 
expected or input Gauss coefficients Pj(t) and recall the diagonal 
matrix R with elements R^j = [n(i) +1) ]5jj : 

g = R p - * (A5) 

Recall that V = R“^ER”^ represents the error covariance matrix for the 
Gauss coefficients which, to the extent that the modeler's expectations 
are realized, measures the (square correlated) uncertainties in the 
Gauss coefficients. Then, in matrix notation, the generalized square 
weighted signal in the radial field coefficients 

[S ( t ) ] 2 = g T Qg = p T R T fiRp = p T R T (RVR T ) _1 Rp = p T V _1 p (A6) 

is the weighted signal in the Gauss coefficients.* 

Now, let Sg^ = g^t) - gi(t 0 ), let Sp.^ s p^t) - Pi(t c ), and omit 
the uncertainty in the initial conditions at time t Q . Relative to t Q , 
the generalized weighted signal in the secular change at t is then 

[AS(t;t 0 )J 2 = 5g T Q5g = 6p T V -1 8p - /(t) [S r *(a,t;t 0 ) ] 2 (A7) 

With initial condition y(t 0 ) = g(t Q ), hence predicted Gauss coefficients 
g{t 0 ) = p(t Q ), the instantaneous generalized square weighted residual 
in the secular change at time t relative to t Q is 

/( t) [8 r *(a,t;t 0 )] 2 = (g - y) T £1 ( g - y) = (p - p) T V' 1 (p - p) (A8) 

The time integral of (A8) from t Q to t f is (A2)-the total generalized 
square weighted residual in the secular change accumulated during the 
interval [t Q ,t^]. 

By replacing W^j/47C with one transforms the objective from 

an attempt to fit the evolution of the weighted radial magnetic flux 
density called for by a time series of geomagnetic field models into an 
attempt to fit the evolution of the scalar geomagnetic potential called 
for by a time series of weighted Gauss coefficients. The former may 
seem more sensible because field components are observable, unlike the 
scalar potential, and because horizontal components of the induction 
equation are not used (0); however, the latter has appreciable merit. 

For example, if the expected Gauss coefficients at time t could be 
obtained from a geomagnetic data vector d, a symmetric data weight 
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(A9) 


matrix M, and a normal equations matrix A such that 
p = (A T MA) _1 A T Md = VA T Md 


then by (A6) 

S(t) 2 = g T tlg = P T V _1 P = d T M T AW” 1 VA T Md = d T M T AVA T Md (A10) 


Now AVA t = A(A t MA) _1 A t = Q implies QMA = A; but Q is not necessarily M" 
becauie A is typically rectangular an3 not invertible. However, if Q 
were M _1 7 then S 2 would equal d T Md-the weighted signal in the data from 
which~the geomagnetic field model could be derived; then (A8) would 
equal the instantaneous signal in the secular change indicated by the 
weighted data, and its time integral would equal the total signal in the 
secular change called for by the weighted data accumulated during the 
interval . 

Unfortunately such data might not exist due, for example, to the 
use of damping, prior bias in deriving the field model, averaging of 
coefficients, or roundoff errors. If such data do exist, then their 
weighting might be suspect; moreover, the data may well contain 
contributions from fields other than the broad-scale core field of 
interest (e.g., crustal fields). The latter seem problematic due to the 
time-varying spatial distribution of geomagnetic survey data. Inclusion 
of such extraneous fields in the objective function seems inappropriate, 
so a suitable modification of Q might be considered. Fortunately, great 
compensation for crustal and external fields can be achieved using the 
correlated data weight matrix technique developed by LANGEL, ESTES & 
SABAKA (1988a, 1989). This kind of approach was used to derive the DGRF 
models for epochs 1945 through 1960 (LANGEL et al., 1988b). Yet it is 
still not clear that general weights are preferable to radial field 
weights-particularly as the correspondence between the weighted signal 
in the data and the weighted signal in the field model is not rigorous. 

As pointed out by G. BACKUS (1987, personal communication), it is 
useful to introduce the quantity 

x k s “k i e i (A11) 

with covariance matrix elements 


Eofx^x^} = Eofo^^e^ejM^ j } = CD^Eofe^ej j 

= ®ki E ij to lj = 8 kl • (A12) 

For a tenth-degree DGRF model at epoch t n the trace of (A12) is 

Tr (Eo{Xj c x 1 ) ) = Eo{Tr (x k x^) } = Eo{x k x k ) = 120 = /(t n ) 

= Eo{ (to k ^e^) © k jej } = Eo{e i E _1 i jej) = Eoje^ijej} (A13) 
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To the extent that the fitting residuals [g^ - Y^] resemble the expected 
uncertainty, the suitably normalized generalized square weighted 
residual at time t n 

[8 r *(a,t n ;t 0 )] 2 = l gi (t n ) - - 7j (t n )]/120 <A14) 

is expected to be unity. This quantity is readily evaluated provided 
ft = ET 1 can be obtained (see Part IC) . 

~0f course, other kinds of geomagnetic field models, notably the 
harmonic spline models of SHURE, PARKER, Sc BACKUS (1982), may not have 
truncated spherical harmonic representations. In such cases, a total 
change of basis functions seems less computationally burdensome than 
transforming the finite dimensional expected error covariance for the 
modeled parameters (say, the harmonic spline coefficients) into an 
equivalent, apparently infinite dimensional, expected error covariance 
for spherical harmonic coefficients of the radial field. 
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ERRATA: PART IA. 

Page 9, line 30. The reference should be upper case and in 
parentheses rather than square brackets. 

Page 22, lines 28-34. The passage "On Earth's surface ... the main 
field. Nevertheless, ..." should be replaced with "Near Earth's surface 
(and apparently within the Earth) the high frequency electromagnetic 
oscillations of solar and terrestrial origin have far greater energy 
density than the main geomagnetic field, so Ampere's law is broken. 
However, ‘ 

Page 25, line 11-12. Note that non-subscripted - E i ^(^t) 

can be a sum of matrices. 
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